Evolving spherical Boson Stars on a 3D cartesian grid 
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A code to evolve boson stars in 3D is presented as the starting point for the evolution of scalar 
field systems with arbitrary symmetries. It was possible to reproduce the known results related 
to perturbations discovered with ID numerical codes in the past, which include evolution of stable 
and unstable equilibrium configurations. In addition, the apparent and event horizons masses of 
a collapsing boson star are shown for the first time. The present code is expected to be useful at 
evolving possible sources of gravitational waves related to scalar field objects and to handle toy 
models of systems perturbed with scalar fields in 3D. 
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I. INTRODUCTION 

Boson Stars (BSs) are localized, regular, spherically 
symmetric solutions of Einstein's field equations, whose 
matter content is provided by a complex scalar field 
with mass and self-interaction 0,1113. Recently Boson 
Stars have appeared in many contexts. Traditionally 
the^were proposed to be alternatives of compact objects 

matter this includes the possibility to emulate 



they ■ 

nil 



; they have been proposed as candidates for dark 



galactic supermassive Black Holes 0; they have even 
been considered as sources of gravitational waves ; and 
they have been used in the field of numerical relativity 
to test evolution formulations of Einstein's equations 
and experience with gauge conditions 0, since these 
objects are very tractable in the sense that they have 
no defined surface and their dynamics do not tend to 
form shocks, and therefore sophisticated techniques -like 
shock capturing methods- are not required. Moreover, 
BSs define sequences of equilibrium configurations which 
serve to test a numerical implementation. That is, for 
a spherical equilibrium configuration the field has to 
oscillate with a constant frequency while the geometry of 
the space-time has to remain static. For general reviews 
in BSs see H, I3- In this paper I focus on numerical 
challenges, and set up a starting point to study the 
astrophysical applications mentioned before. Instead 
of developing the original mean field approximation 
to calculate the expectation value of the stress-energy 
tensor cornponents of a quantum scalar field as done 
originally P, 0) the field is simply considered to be 
classical as usually done when studying the evolution of 
these systems. 



For several reasons, such systems have only been 
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evolved for long times (long enough to observe the 
oscillations of a perturbed star) in spherical symmetry 
using ID codes \wL llll ll^|. The main reason for this is 
without doubt the difficulty in formulating the evolution 
equations for the space-time from Einstein's equations 
when using non-spherical symmetry. In spherical sym- 
metry these evolution equations are replaced with the 
slicing condition and the constraints (see for instance 
[Ulia)- The result ing ODEs can be integrated at each 
time step during the evolution, which is computation- 
ally cheap and allows a constrained evolution under 
control for very long periods of time. Other problem in 
higher dimensions is the difficulty to achieve the same 
resolutions used in ID simulations due to the lack of 
computational resources, since the memory requirements 
in 3D unigrids forces one to sacrifice either the size of the 
physical domain or the numerical resolution in regions 
where the gravitational field still strong. Moreover, the 
number of variables in higher dimensions is also bigger 
than in ID. 

The practical concern is that in 3D the constraint 
equations are elliptic, which one does not want to 
solve at every time step due to the extreme hardware 
requirements needed to carry out a significantly long 
simulation. Therefore it is usually preferable to use 
unconstrained gauges and evaluate the constraints just 
to monitor the accuracy and convergence of the solution 
one is calculating. Even though unconstrained evolutions 
of Einstein's equations are still under study, it is known 
that some formulations work better than others when 
using certain gauge conditions. In this paper the BSSN 
formulation of general relativity is considered i which 
has shown to be more stable for longer evolutions than 
the ADM formulation for a number of physical systems 
and presents different stability properties Q (for 
test beds and comparisons under different gauges see 

M)- 

The Lagrangian density from which the Einstein and 
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Klein-Gordon equations describing BSs are derived reads: 



II. NUMERICAL METHODS 
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where R is the Ricci scalar, (/> is a complex scalar field 
minimally coupled to gravity, 0* is its complex conju- 
gate and V the potential of self-interaction of the field. 
The resulting Einstein's equations are G^i/ — 8ttGT^„, 
provided the stress energy tensor is 



T„ 



lid. 



(2) 

The equation satisfied by the scalar field can be obtained 
either from the Bianchi identities or by directly using the 
variational principle on which for the present case 
reduces to the Klein-Gordon equation on the space-time 
background. 



□ - 



dV 



= 



(3) 



where n(j) — 
sponding line element: 
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where latin indices stand for the spatial coordinates 
and X indicates dependence on all spatial (cartesian) 
coordinates, a is is the lapse function, /?' the shift vector 
and 7ij is the 3-metric in the 3-1-1 decomposition of the 
space-time. As I am going to deal only with spherically 
symmetric Boson Stars, a zero shift is used for stable 
equilibrium configurations, and a simple but sophisti- 
cated driver involving a non-zero shift will prove to be 
necessary in order to follow the evolution of unstable 
configurations. The general form of the potential to be 
studied here is 1^(101) = ^to'I^P -I- j\4>\^, which has 
been the most studied case with ID codes ll^ Hit IT^ . 
among a variety of potentials studied 0,0,0] and oth- 
ers including very general non polynomial potentials |19| | . 

With this notation in mind in section^ the numerical 
techniques used in the simulations shown are described. 
In section liiil physical tests are presented and a set 
of production runs from which the frequencies of the 
oscillation modes found in ID for several configurations 
are reproduced with the present 3D code. In section 
IIVI an example of various unstable BSs is shown. A 
trick to excite non-radial modes is presented in section 
Ivl through the manipulation of the refinement of the 
grid in different directions. Section IVII refers to the 
restrictions of the present code. Finally, some comments 
are addressed in section IVIll 



A. Evolution in time 

The scalar field is described by its real and imaginary 
parts, such that 0(x, t) = 0i(x, t)+i(j>2{x, t). I define four 
new real variables based on combinations of derivatives of 

1/2 

the scalar field as follows: tti -|- i7r2 := ^^^-^{dt4> — P^dicj)), 
and ipia + iip2a '■= da4> for a = x,y, z, where a is still the 
lapse function of the space-time and 7 is the determinant 
of the three metric 7^ . In terms of these new variables 
the KG equation (O can be rewritten as a set of first 
order evolution equations: 



dtipia 

dtT^2 



9. («7'/VV^ij+/3Vi 
9. (a7'/'7''V^2, + /3V2) 



(5) 



(7) 



1 1/2 dV ^/c 



where the summation convention is assumed over latin 
indices. Evidently the set H5I8() is a system of evolution 
equations that provides the derivatives of the scalar 
field at each time slice, from which it is possible to 
reconstruct the scalar field itself using the definition 
of TTi and TT2- When solving the KG system above, it 
is necessary to solve the evolution Einstein's equations 
in the BSSN formulation as described in reference [l3| . 
Basically the BSSN formulation uses the evolution vari- 
ables: * = ln(det7y)/12, % = e-^*^,^-, K = -/'^K,j, 

--!*(/«:„ - 7yi^/3) and = 75*Tk, instead of 



Aij — e 

the usual ADM variables 7ij and Kij. The evolution 
equations for the new variables are described in a 
number of references, e.g. [T^ . 

The integration of the Klein-Gordon equation I|5I8|) 
and the BSSN variables is performed using the method 
of lines with second order centered differencing in space. 
For the time integration a third order Runge-Kutta 
(RK3) algorithm is used 20] . 



B. Gauge choice 

In the case of stable equilibrium configurations the evo- 
lution of the gauge was carried out using zero shift and 
the 1 -I- log slicing condition which is a particular 
method according to which the lapse evolves as 



dta = -2a{K - Ko) 



(9) 
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where K is the trace of the extrinsic curvature and 
Kq = K{t = 0) {Kq—Q in the present case). Due to 
the time independence of the metric one knows a priori 
that Kij = —^dtjij — 0, i.e. that K has to be zero 
for equihbrium configurations of boson stars. However 
this is only valid in the continuum limit. In practice, 
the discretized version of the equations will prevent the 
extrinsic curvature from being exactly zero. Nevertheless 
it was observed that the components of the extrinsic 
curvature converged to such value. 

In order to evolve unstable equilibrium configurations 
the 1 + log slicing condition was also used, but it was nec- 
essary an additional hyperbolic gamma-driver shift con- 
dition This type of conditions are necessary after 
a black hole has formed and prevent the grid stretching 
effect followed by an unphysical growth of the horizon. 
As described in |22| . such gamma-drivers are diverse, but 
for the spherically symmetric case I deal with here, the 
following hyperbolic condition was used: 

dtP' = ^aPB' , dtB' = dtV - TjB' (10) 

where p and 77 are coefficients controlling the gauge 
speeds and certain damping added to the shift respec- 
tively. 

C. Exterior Boundary Conditions 

A modified version of radiative boundary conditions is 
applied to the BSSN variables and to the defined first 
order variables 0, ip and tt. For all these variables it is 
assumed that they behave as a constant plus an outgoing 
radial wave at the boundaries, plus an extra term that 
helps to control initial transient effects, that is 

f{r,t) = fo + u(r-tyr + h{t)/r^ (11) 

where r = |r|, /o is one for the lapse and the diagonal 
components of the three metric and zero otherwise, h{t) 
is a function of t and n an unknown power. The reason 
to choose the outgoing radial wave boundary condition 
is that it has proved to be very efficient at absorbing 
waves. The differential equation in Cartesian coordinates 
corresponding to the expression above is 

^aj + 9./-H^(/-/o)^^. (12) 

For a given n, the unknown function h'(t) can be 
evaluated at a point away from the boundary, and 
substituted in the previous equations to evaluate 
any dynamical variable at the boundary p^ . It was 
found that when the potential of the scalar field has 
zero self-interaction (A = 0), then switching off the 



extra term {h = 0), that is, using standard radiative 
boundary conditions, the evolution behaves properly. 
Nevertheless, when the self-interaction is non zero, 
it was found that a bigger power (n = 3 or n = 5) 
prevents a pulse on the field variables coming from 
the boundaries. It is worth mentioning that instead 
of using a sponge as used in spherically symmetric 
evolutions with A 7^ [TH . the modifications mentioned 
here sufficed to evolve the shown systems with second 
order accuracy for a number of crossing times. Finally 
it is important to stress that in the case of spherical 
symmetry these BCs are appropriate. More care is re- 
quired when systems with other symmetries are studied, 
since this approximation is only valid if one considers 
that the front signals propagating outwards are spherical. 



D. Interior Boundary Conditions 

After an apparent horizon has been found for the first 
time, a lego sphere (the set of points of a cubic grid 
contained inside a 2-sphere) is excised from its interior. 
Excising is a common practice when evolving black hole 
space-times 21, 2^ [2^ |^ and has shown to be an 
essential strategy when evolving matter fields as shown 
recently in psl [2q| ; the main aim of using excision is that 
the code is not forced to calculate quantities in a region 
close to a singularity, where gradients are expected to be 
inaccurate. In the present case one expects the matter 
field to fall into the resulting horizon in a finite time, 
thus the excision region is allowed to resize according 
to the size and shape of the apparent horizon. For this 
purpose the apparent horizon finder in was used 
every time step. A buffer zone of five points between the 
excision region and the apparent horizon is kept. 

The interior boundary is the set of points delimiting 
the excision region. As apparent horizons lie inside event 
horizons, all the signals propagate towards the excision 
region and so do the errors. This is the reason why the 
time derivative of all the evolved variables is set to be 
zero at the the inner boundary. The implementation con- 
sists of calculating the time derivatives of the variables 
one point away from the excision region in the direc- 
tion determined by the normal; this time derivative is 
then used to evolve the variable on the excision boundary. 



III. TESTS ON THE STABLE BRANCH 

A. Initial Data 

In order to test the evolution code, the first important 
issue is setting up initial data that work in ID. For this, 
the well studied case of spherically symmetric equilib- 
rium configurations is chosen, where the scalar field os- 



cillates with constant frequency whereas the geometry in 
the chosen gauge remains time-independent. Einstein's 
equations are solved on a spherically symmetric back- 
ground by assuming that the space time is static and that 
the scalar field has the form (j){r,t) = (/3(r)e"^*, so that 
all the components of the stress-energy tensor are time- 
independent. Assuming ds^ = —a{r)'^dt'^ + a{r)^dr^ + 
r^dil^ one arrives at the equations 
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which is a set of coupled ordinary differential equations 
to be solved under the conditions a(0) — 1, (p{0) finite 
and (pf{Q) = in order to guarantee regularity and 
spatial flatness at the origin, and (p{oo) = ipf(oo) = 
in order to ensure asymptotic flatness at infinity as 
described in H E [ill |2|. Here k = SttG. The 
solution was calculated using the standard rescaled 
variables ip — V^nGip, r — mr, t — cut, a — —a and 



A 



so that the coordinate time is given in units 



A 

of l/cj and the distance in units of 1/m. In Fig. Qthe 
sequence of equilibrium configurations corresponding to 
different self-interaction values is shown, and will serve 
to illustrate the functionality of this 3D code. 

Once the system (|13f) has been solved with sufficiently 
high resolution, the straightforward identification of 
the ID variables with the 3D varia bles at in itial time 
is as follows: a{x,0) = a{r), ^nVrlxTO) — a{r), 
(/)i(x,0) = (fir), 02(x,O) = 0, 7ri(x,0) = and 
7r2(x, 0) = ^fir)- This can be easily understood just by 
recalling the ansatz (j){r, t) — (/?(r)e** in rescaled units. 
Then the 3D cartesian grid is populated by copying the 
values of such variables to the nearest neighbor grid 
points of the 3D domain, and applying the coordinate 
transformation from spherical to cartesian coordinates. 
Since it is known a priori that the space time should be 
time-independent I set the components of the extrinsic 
curvature to zero initially. Then the unconstrained 
evolution is started and the geometry and the field 
evolve according to their corresponding equations. 



B. Convergence 

The convergence of a non vacuum system is a major 
issue for many reasons. The usually applied boundary 
conditions to the geometrical quantities correspond to 
outgoing waves and moreover the boundary conditions 
for the scalar field with a potential term have no 
reported solution up to now for a 3D grid, so that 




(p(0) 



FIG. 1: The sequences of equilibrium configurations are 
shown for different values of the self-interaction coefficient A. 
The configurations to the left of the respective peaks are sta- 
ble configurations, and those to the right are unstable, which 
means that either they could collapse to form black holes or 
disperse away under infinitesimal perturbations, depending 
on the sign of the binding energy (see below). The triangles 
down, the stars and the circles correspond to the configura- 
tions evolved to trace the radial modes of oscillation of stable 
configurations. The filled squares represent four particular 
cases studied in depth to show the capabilities of the numer- 
ical techniques. 



the boundary conditions are not necessarily consistent 
with the evolution equations; these difficulties are 
not particular to the present problem, but are also a 
concern when solving the vacuum Einstein's equations 
or relativistic hydrodynamical fields, because these 
conditions could allow constraint violating modes to 
come into the domain and destroy the simulations. 
Nevertheless it was found that the modified Sommerfeld 
boundary conditions mentioned in the previous section, 
work well for the variables describing the field and the 
geometric quantities at once for quite a long time. In 
Fig. 121 it is shown the convergence of the L2 norm of 
the violation of the Hamiltonian constraint produced for 
two initial configurations. Second order convergence is 
manifest indicating the accurate and stable behavior of 
the simulations for dozens of crossing times. 

One intriguing issue has to do with the time inde- 
pendence of the metric, which actually oscillates around 
or close to a certain fixed shape. In fact, due to the 
considerable amplitude of such oscillations it is possible 
to imagine that probably the evolved configuration 
corresponds to another configuration with different 
properties than those of the initial data. In Fig. Otwo 
different cases that are representative are shown: in 
the first case, the maximum of the metric function ^xx 
oscillates around the right value calculated with the 
ID code, which is assumed to be the correct value; in 
the second case the maximum of the metric component 
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A=30, 0(O)=O.O6, M=0.8225 (M/mJ 
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FIG. 2: Convergence of the relative L2 norm of the hamilto- 
nian constraint for two configurations; the normalization fac- 
tor is the T't component of the stress energy tensor. The 
resolutions used were Axyx — 0.675 for the 40^ box and 
Axyz = 0.3375 for the 80^ box in each of the configurations. 



oscillates, but not around the expected value, instead 
with smaller values all the way. The fact is that in both 
cases the oscillating values of max{'jxx) converge to a 
constant value that coincides with the one calculated 
at initial time. The second order convergence of this 
quantity indicates that the desired configuration is the 
one being evolved. 



C. Radial Modes 

With high spatial resolution it is quite easy to obtain 
very small variations of the metric functions in time, but 
since one is dealing with a coarse 3D grid the situation is 
slightly different. Instead, due to the errors introduced 
when populating the 3D grid with the initial data and 
those due to the discretization, one is dealing with 



perturbed configurations and thus it is possible to see 
only how the geometry converges to a fixed one when the 
spatial resolution is increased. Nevertheless I take ad- 
vantage of this fact now to study oscillating radial modes 
as done in ^Jli but in the present case using the 3D code. 

When studying the ID problem with very refined 
grids it was necessary to set up an explicit pertur- 
bation by adding particles or kinetic energy to the 
initial equilibrium configuration in order to study the 
oscillations of the space time ^- Taking advantage 
of having a coarse grid, one gets for free the errors 
introduced when setting up the initial data, plus the 
discretization errors due to the numerical methods, act 
as a spherical perturbation, since an equally spaced grid 
is used in all directions, so that the evolution starts 
with perturbed configurations that converge to a strict 
equilibrium configuration, and it is possible to follow 
the oscillations of the metric functions. An indication 
that the perturbations applied are spherical consists in 
showing the amplitude of the Zerilli function extracted 
from the evolutions for the case of an unequally spaced 
grid in the z-direction as shown below. 

In Fig. ^ appear the results found for equilibrium 
configurations evolved for three different values of 
self-interaction. Each point represents a run on a 40'^ 
(the coarse grid here) in octant mode of one of the 
equilibrium configurations indicated in Fig^ which was 
run up to i > 1000 in coordinate time and stopped 
because of two factors: a) even if the evolution algorithm 
used is less dissipative than the usual second order ones, 
dissipation starts being noticed and b) the phase of the 
scalar field also starts to be shifted and begins to be less 
good. The frequency of the modes shown was calculated 
as usual / = l/T/a{oOgrid) where T is the period of 
the maximum of the metric functions and oogrid is 
the outermost point of the grid. The curves shown in 
Fig. ^ mimic those found for these same configurations 
obtained with a ID code important difference 

between the method used in and the used here 
is that in the former case the system was perturbed 
with a shell of particles falling towards the star, and 
in the present case the discretization error is used as 
the perturbation. In the first situation the metric 
functions oscillate with an amplitude that approaches 
exponentially towards their values corresponding to the 
unperturbed configuration, and in the present case no 
extra particles are added, which implies that there is no 
extra matter to be ejected, and no exponential decay of 
the amplitude of the oscillations is expected, instead, 
the discretization error is present over all the numerical 



Nevertheless this author has measured the same effects without 
the need of perturbing the number of particles and just taking 
care of the perturbations due to the discretization errors in a 
personal ID code. 
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FIG. 3: On the left, the j^x metric function is shown for different systems evolved on an 80'^ grid and the highest resolution 
mentioned in the previous figure. On the right is shown the shape of the maximum of y^x in time and its second order 
convergence to the constants max{'yxx) = 1.14276 and max{'yxx) ~ 1.12525, which correspond to the values for each of the two 
configurations calculated at initial time with high resolution using a ID code. 



domain at every time. 



IV. THE UNSTABLE BRANCH 

Unstable boson star configurations are equilibrium 
configurations characterized by the fact that even using 
very fine tuned initial data and very high resolution, 
the finite differencing error in the calculations suffices 
to make the star collapse or disperse. In order to 
analyze the physical properties of unstable BSs I study 
in particular the case A = 0. An important feature of 
the unstable branch of BSs is the threshold related to 
the sign of the binding energy Eb — M — Nm, where 
N is the number of particles N = J and M is the 
mass measured by an observer at infinity, li Eb < the 
system is expected to collapse and form a black hole, 
and otherwise the system is unstable due to the excess 
of energy translated into kinetic energy, and (like in a 
fission process) it should be dispersed into free particles 
at infinity 28] , and this is precisely what is confirmed in 
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FIG. 4: The frequencies of oscillation of miLK{'yxx) for several 
equilibrium configurations is shown. It was found to be in very 
good agreement with those plots found with a ID code for the 
same systems and others in reference Jjj . Here I just present 
some cases corresponding to certain values of A = X/AnGm^ . 
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these simulations^. The hmiting situation Eb = occurs 
for the case A = when 0(0) = 0.54. This threshold 
splits the unstable branch of equilibrium configurations 
into two regions: one corresponding to configurations 
with central field value 0.27f < (j){0) < 0.54, i.e. 
between the critical point (/)(G) = 0.271 separating the 
stable and the unstable branches and the threshold 
of zero binding energy (j>{Q) = 0.54; the fate of the 
configurations belonging to this region is the formation 
of a black hole. The second region corresponds to con- 
figurations with (j){0) > 0.54, which should disperse away. 

In Fig. |S1 a typical density evolution is presented 
for the case of a collapsing BS, both before and after 
the Black Hole is formed. After finding an apparent 
horizon for the first time, the excision is applied to all 
the variables that are being evolved and those serving 
to monitor the system (e.g. the energy-density). In the 
same figure the mass of the apparent and event horizons 
of the resulting black hole are shown. The discrepancy 
between the mass of the apparent and event horizons 
and the original mass of the system is within 0.4% after 
most of the matter has been absorbed, around t ^ 90. 
The event horizon is a concept that depends on the 
global structure of the space-time, therefore it can only 
be tracked in a post-process step after the simulation was 
finished. The event horizon finder j29| works by evolving 
null surfaces backwards in time using the information 
stored during the space-time evolution. Two initial 
surface guesses are necessary to start the event horizon 
tracking. One guess is chosen in a region outside and 
other inside the expected location of the horizon. These 
guesses were chosen based on the apparent horizon 
radius. The event horizon is at the end the surface the 
two initial guesses converge to. This is the reason why in 
Fig. [SI two different event horizon masses are shown. It 
can be seen that by the end of the simulation (the start 
of the event horizon tracking) the mass of the two initial 
surfaces does not coincide, and that around t — 105 the 
two measurements converge to a single one: the event 
horizon. 



Based on these two important tools it is shown in 
Fig.Elthat the ratio of polar to equatorial circumferences 
is one over the evolution, indicating that the horizons 
are spherical, despite the lego-type shape of the excision 
boundary. This ratio will be fundamental in determining 
the physical parameters of the resulting black hole in the 
case of evolutions with fewer symmetries. E.g. from the 
oscillations of Cr = Cp/Ce it is possible to obtain the 
angular momentum of a black hole resulting from the 
collapse of very general configurations (see e.g. |30ll3l|'). 
Finally, in Fig. [7| the evolution of for a dispersing BS 
is presented, where the ejection of matter is manifest. 



V. EXCITING NON-RADIAL MODES 

Since the systems presented here have spherical 
symmetry, no gravitational wave signals are expected. 
Nevertheless as an exercise and in order to show that 
the modes studied in section HTll were actually radial, I 
take advantage of the freedom one has in 3D to apply a 
simple trick just by choosing different spatial resolutions 
Ax = Ay ^ Az with the consequent effect that the 
system is being perturbed in the z direction, and thus 
gravitational radiation is obtained in the fundamental 
mode simply due to discretization error. The pertur- 
bation consisted in using a grid size 40 x 40 x 35 in 
the x,y,z directions respectively and resolution such 
that Ax = Ay, Az was chosen so that the range in the 
three directions was the same. In Fig. [S] the waveform 
for the fundamental mode I = 2, m = are shown 
and compared with the results found for an isotropic 
grid. The wave forms were extracted with the methods 
described in [s^]- The convergence of the wave forms 
is second order, and as expected it converges to zero, 
since in the continuum limit it is expected the grid to 
be isotropic and no signal in this mode could come out 
from a spherically symmetric source. In this way it is 
also shown that the modes analyzed in the section IIIII 
were actually radial. 



^ In (ir| it was found that independently of the sign of the bind- 
ing energy, the configurations belonging to the unstable branch 
coll apse into black holes for the case A = studied here. In 
llll however, it is confirmed that for Eb > in the case of ex- 
cited boson stars the configurations are dispersed away. In the 
present 3D case, due to the limitations on the numerical domain, 
it was not possible to confirm is after an expansion (far away from 
the 3D box) the system should have recoUapsed to form a black 
hole, or if a migration to a very dilute configuration ocurred. In 
any case, it is a coincidence that stars with (/"(O) > 0.54 simply 
dispersed away, whereas those with 0.271 < <l>{0) < 0.54 simply 
collapsed without any noticeable ejection of matter. An essential 
difference to be kept in mind is the nature of the perturbations 
applied here, where the number of particles is not being modified 
at all, and the perturbation is applied to the system across the 
whole numerical domain through discretization error 



VI. FEATURES OF THE CODE 

A. Restrictions 

The period of time that the convergence is preserved 
acts as an indicator of the consistency properties of the 
scheme, formulation, and boundary conditions working 
all together. The runs were stopped at the times 
shown in the plots because: a) the convergence was not 
manifestly second order with the passage of time, b) 
the period of the scalar field drifted after a finite time. 
These problems could be sorted out just by using higher 
resolution in time; all the simulations were carried out 
using a CFL factor of At /Ax = 0.05 (which means 
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FIG. 5: In the top panel the evolution of the energy density is shown. At the first stage of the evolution (between t = and 
20) the density remains time independent; nevertheless after t ~ 40 the distribution starts concentrating close to the origin, 
and after f ~ 64 the system collapses and an event horizon is formed. Then around t ^ 68 an apparent horizon is found, after 
which excision is applied. At later times the remaining matter outside the horizon keeps falling into it. In the bottom panel 
the mass of the apparent horizon (Mah), the event horizon (Meh) and the initial total mass (Minit) are compared. After 
t = 90 and towards the end of the simulation the discrepancy is less than 0.4% between Meh and Mi„it; a zoom in would 
show that the AIah is always less than Meh- The point marked with a cross indicates the moment at which an apparent 
horizon appears, after which, the increase in its mass is due to the accretion of the scalar field remaining outside. The two 
lines indicating the mass of an event horizon should be understood as follows: an interior and an exterior initial guesses for 
tracking the event horizon are chosen, thus it is expected they converge to two different null surfaces at the beginning of the 
tracking (that is, at the end of the simulation); however we know that when these null surfaces coincide the event horizon has 
been hunted; therefore the event horizon starts being tracked at around t = 105. The configuration corresponds to A = and 
(t>{0) ~ 0.4, which is a typical case of the simulations collapsing into a black hole. The gauge parameters in IIOII are p — 2 and 
■q — 0.25. The evolution was performed with resolution Axyz — 0.07875 in an 80"^ box in octant mode. The simulation was 
stopped at around £ ~ 117 due to an instability coming from the inner boundary. 
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for the different horizons, and it is manifest that their shape 
is spherical within an error of 0.2% towards the end of the 
simulation. 
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FIG. 7: The energy density profile of the configuration A = 
and (ji{0) — 0.6 is shown at difi'erent stages. Due to the re- 
stricted domain it could happen either that the system should 
recoUapse after traveling a distance outside the box to form 
a black hole as claimed in |T(i^ . or it is confirmed the fact 
that for configurations with positive binding energy the sys- 
tem truly disperses away and never recoUapses. This contro- 
versy is expected to be solved once the present code is able 
to use notably bigger domains with the help of mesh refine- 
ment techniques. The configuration corresponds to A = 
and 0(0) = 0.6. The gauge parameters and resolution are the 
same as those used for the collapsing system shown above. 



that the runs required tens of thousands of iterations), 
but convergence and period could be preserved just 
by decreasing such number. Another reason to choose 
such a small value of the CFL factor was that it helps 
at diminishing dissipation after a very long time of 
evolution. 

Another issue is the problem of resolution. During 
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FIG. 8: The wave form and its convergence is shown for the 
fundamental mode with a detector located at a; = 30(l/m). 
The perturbed system was A = 0, cl>iiO,0) = 0.09 and M = 
0. 51533 {M pi /m). The solid line around zero is the amplitude 
of the wave form when the system is evolved on the isotropic 
grid. 



the evolutions presented in this manuscript there were 
at most a couple of dozens of points (using the low 
resolution) covering the volume containing 95% of the 
total mass of the star, and twice as much outside such 
region. In consequence the gravitational wave detectors 
had to be placed rather close to the BS. But even so, the 
results are good. The ideal situation for more general 
configurations would be: high resolution in the central 
part of the grid and a physical domain big enough to 
allow one to place detectors far away from the source. 
Thus practically means that the detector would be 
located in a Schwarzschild-like region thus providing 
cleaner signals. 

Using a uniform grid works fine for spherical sym- 
metry, but presents a dilemma in other cases: use high 
resolution to cover the central part and then forget 
about obtaining high quality waveforms, or conversely, 
sacrifice the resolution at the center and pull the grid 
far away in order to have a clean gravitational wave 
signal, but produced by a not accurately evolved source. 
The solution to this problem would be provided by the 
implementation of mesh refinement [33 . l3^ l35| . and 
will be useful when dealing with systems under more 
interesting symmetries. 



B. Restrictions on the physics 

The phase of the field is an important indicator of the 
control one has on the simulation, because if the field 
is out of phase or the period of oscillation of the field 
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are not the right ones, then the effects will be noticed 
in the geometric quantities. Among the experiments 
carried out, it was observed that when using a three step 
iterative Crank-Nicholson (3ICN) [s^ or second order 
Runge-Kutta (RK2) methods to evolve the system, the 
period of the scalar field changed after a while {t ~ 200) , 
whereas the third order Runge-Kutta preserved the 
phase and the period, although a modified version of 
ICN [37] proved to work properly as well. 



The dissipation is another issue one has to deal with 
when evolving matter. By observing the amplitude 
of the scalar field it was found that after t ~ 1000 
the 3ICN provided an amplitude of the field ~ 10% 
less than the original one, while with the RK2 the 
amplitude of the field it was 15% bigger than the right 
value for the resolutions used here, which effects were 
evident in the measurement of physical quantities. The 
amplitude of the field after t ^ 1000 was less than the 
original by 0.3% using the coarse resolution. Thus, the 
RK3 was the one providing the most confident results 
in terms of convergence for longer time, dissipation, 
phase and period of oscillation of the solutions. It was 
also essayed the fourth order Runge-Kutta scheme, 
however it required more time iterations and no major 
improvement in accuracy was found. RK3 showed to be 
optimal in performance and accuracy for the purposes 
of the runs shown here. 



In the case of stable configurations, the duration 
time scale of the simulations was enough to observe 
the oscillations of a wide range of stable configurations. 
However, the expected period of oscillation for very 
diluted configurations (those to the left in Fig. 0} is 
bigger if the mass is smaller. Therefore, it was not 
possible to study those configurations within the ranges 
of evolution times achieved with the present resolution 
and evolution times at the moment. 



In the case of unstable configurations it is not expected 
to observe a constant amplitude of the scalar field, in- 
stead the squeezing of the scalar field, which is confined 
to a small region, implies an increase in the amplitude 
of the field and energy density whereas it is more con- 
centrated in a small region where nodes can be observed 
before the formation of an apparent horizon. After exci- 
sion is applied, the gradients of the scalar field become 
important in the region close to the inner boundary, and 
some unstable modes develop at some point. We are 
certain these modes are not numerical artifacts because 
increasing the resolution they last more time to show 
up. I expect to report about this point in further reports. 



VII. COMMENTS AND CONCLUSIONS 

I have presented a 3D code that is able to keep 
under control the unconstrained evolution of spherically 
symmetric Boson Stars, both stable and unstable using 
the BSSN formulation of General Relativity, the 1 -I- log 
slicing condition, a coarse grid and simple discretizations 
on a uniform grid and well known numerical methods to 
evolve in time. 

Stable configurations were evolved during several 
crossing times with second order convergence of the 
hamiltonian constraint and some characteristic prop- 
erties of the geometry, like the convergence of the 
maximum of the radial metric function to the ap- 
propriate value. The process of collapse of unstable 
configurations was followed for quite a long time of 
around t ~ 80M, after the BH was formed, with error 
in the mass of the horizons of less than 0.4% in several 
cases studied. The configurations that dispersed where 
followed for a very long time and no surprises were found 
when evolving the resulting flat space-time. 

It has been shown that these evolution techniques 
work properly when the IVP is not solved in the grid 
used for the evolution. This procedure will permit 
one to distinguish future problems intrinsically related 
to the solution of the initial value problem directly 
on the 3D grid js^. I also decided to use simple 
uniform grids in order to set up test beds for future 
Boson Star simulations without involving coordinate 
transformations |^ that will be used for further re- 
search or mesh refinement that will be helpful as well [33 ■ 

It is possible to conclude by stating that it has been 
possible to reproduce in a 3D cartesian grid, the basic 
results found with ID codes for spherically symmetric 
Boson Stars, which are related to infinitesimal pertur- 
bations due to the discretization error. In addition, 
taking advantage of the freedom one has in 3D grids to 
choose the discretization error in different directions, it 
was possible to extract a wave form that illustrates the 
type of signals one could expect from BSs and showed 
that the signal obtained when exciting the fundamental 
mode, presents a quick damping as predicted for axial 
perturbations of Boson Stars 40| . 

A comment about the late development of a 3D code 
for boson stars is in turn. Under the ADM formulation 
of General Relativity it was rather difficult to evolve 
a Boson Star in 3D for the periods of time shown here 
(see e.g. 0)- The main reason preventing such long 
evolution was the ADM formulation itself, which showed 
to be inefficient at evolving a garden variety of physical 
systems On the other hand, at the present time one 
accounts with several new gauge conditions allowing the 
control of the simulations. A few years ago considerable 
effort was invested to evolve Boson Stars in 3D that did 
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not succeed only because of the formulation used |42ll43| . 
although previous studies using new formulations were 
also explored It is necessary to point out that now 
it is possible to study in depth the phenomenology of 
Boson Stars in full 3D with the new techniques and tools 
that are emerging in the field of numerical relativity. 
Boson Stars could still be systems helpful at testing new 
numerical techniques, but most important they can be 
considered objects whose existence could eventually be 
probed if they are thought to be sources of gravitational 
waves. 

This code was written as a thorn in the Cactus frame, 
and can be plugged into the variety of analysis tools and 
drivers available in the Cactus toolkit l4a . 
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